%% Comparative Statics
% Compute and plot comparative statics results

myCluster = parcluster('local');
myCluster.NumWorkers = 24; 
saveProfile(myCluster);  
parpool(24)

%% Timestamp

timestamp = now;

%% DIRECTORY


%% State Space
n_grid= 51;

phi_min = 0;
phi_max = 1;
PHI = linspace(phi_min,phi_max,n_grid)';
dphi = (phi_max-phi_min)/(n_grid-1);

y_min = -1.5;
y_max = 0;
Y = linspace(y_min,y_max,n_grid)';
dy = (y_max-y_min)/(n_grid-1);

%Initialize Matrices
PHI_mat = ones(n_grid,1)*PHI';
Y_mat = Y*ones(1,n_grid);
 

%% Grid params
grid_params = {n_grid, Y_mat, PHI_mat, dy, dphi};


%Numerical parameters
max_iter = 1000;
D = 10000000;
thr = 10^(-9);

num_params = {max_iter, D, thr};

%Control space
M = 5;
eta_min = 0;


%% Base Paramters
top_ret_b = 0.02; 

% Preference Parameters
rho_b = 1/3; 
r_b = 0.02;
delta_b = 0.05;
ra_b = delta_b/(1-rho_b);

% Return Parameters
alpha_b = 0;
sigma_b = 0.18;

% Moral Hazard parameters
lambda_b = 0.95;

%Parameters for Mature Phase
eta_max_b =  top_ret_b/sigma_b;



% Check condition
eta_max_b*sigma_b < lambda_b*sigma_b*rho_b*sqrt(2*(1-rho_b))*sqrt( (ra_b - r_b)*(1-rho_b)/(rho_b)  )

% Parametes to feed to external functions

params_b = [rho_b, r_b, ra_b, alpha_b, sigma_b, lambda_b, M, eta_min, eta_max_b];

%% Benchmark
[~, C_b, Beta_b, K_b, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, vol_c_b, ~] = ...
    s_benchmark_system_solver(params_b, grid_params, num_params, timestamp);

epsilon_C_b = (vol_c_b./C_b + Beta_b)/sigma_b;

%% Eta
top_ret_h = 0.021;
eta_max_h = top_ret_h/sigma_b;

top_ret_l = 0.019;
eta_max_l = top_ret_l/sigma_b;


% Check condition and solve
eta_max_h*sigma_b < lambda_b*sigma_b*rho_b*sqrt(2*(1-rho_b))*sqrt( (ra_b - r_b)*(1-rho_b)/(rho_b)  )
params_eh = [rho_b, r_b, ra_b, alpha_b, sigma_b, lambda_b, M, eta_min, eta_max_h];
[~, C_eh, Beta_eh, K_eh, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, vol_c_eh, ~] = ...
    s_benchmark_system_solver(params_eh, grid_params, num_params, timestamp);

epsilon_C_eh = (vol_c_eh./C_eh + Beta_eh)/sigma_b;


dC_dE = (C_eh - C_b)/(eta_max_h - eta_max_b);
dB_dE = (Beta_eh - Beta_b)/(eta_max_h - eta_max_b);
dK_dE = (K_eh - K_b)/(eta_max_h - eta_max_b);
dEps_dE = (epsilon_C_eh - epsilon_C_b)/(eta_max_h - eta_max_b);

%% Lambda
lambda_h = 0.96;

% Check condition and solve
eta_max_b*sigma_b < lambda_h*sigma_b*rho_b*sqrt(2*(1-rho_b))*sqrt( (ra_b - r_b)*(1-rho_b)/(rho_b)  )
params_lh = [rho_b, r_b, ra_b, alpha_b, sigma_b, lambda_h, M, eta_min, eta_max_b];
[~, C_lh, Beta_lh, K_lh, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, vol_c_lh, ~] = ...
    s_benchmark_system_solver(params_lh, grid_params, num_params, timestamp);

epsilon_C_lh = (vol_c_lh./C_lh + Beta_lh)/sigma_b;

dC_dL = (C_lh - C_b)/(lambda_h - lambda_b);
dB_dL = (Beta_lh - Beta_b)/(lambda_h - lambda_b);
dK_dL = (K_lh - K_b)/(lambda_h - lambda_b);
dEps_dL = (epsilon_C_lh - epsilon_C_b)/(lambda_h - lambda_b);

%% r
r_h = 0.01;

% Check condition and solve
eta_max_b*sigma_b < lambda_b*sigma_b*rho_b*sqrt(2*(1-rho_b))*sqrt( (ra_b - r_h)*(1-rho_b)/(rho_b)  )
params_rh = [rho_b, r_h, ra_b, alpha_b, sigma_b, lambda_b, M, eta_min, eta_max_b];
[~, C_rh, Beta_rh, K_rh, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, vol_c_rh, ~] = ...
    s_benchmark_system_solver(params_rh, grid_params, num_params, timestamp);

epsilon_C_rh = (vol_c_rh./C_rh + Beta_rh)/sigma_b;

dC_dR = (C_rh - C_b)/(r_h - r_b);
dB_dR = (Beta_rh - Beta_b)/(r_h - r_b);
dK_dR = (K_rh - K_b)/(r_h - r_b);
dEps_dR = (epsilon_C_rh - epsilon_C_b)/(r_h - r_b);



%% rho
rho_h = rho + 0.01;
ra_h = delta_b/(1-rho_h); 

% Check condition and solve
eta_max_b*sigma_b < lambda_b*sigma_b*rho_h*sqrt(2*(1-rho_h))*sqrt( (ra_b - r_b)*(1-rho_h)/(rho_h)  )
params_rhoh = [rho_h, r_b, ra_h, alpha_b, sigma_b, lambda_b, M, eta_min, eta_max_b];
[~, C_rhoh, Beta_rhoh, K_rhoh, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, vol_c_rhoh, ~] = ...
    s_benchmark_system_solver(params_rhoh, grid_params, num_params, timestamp);

epsilon_C_rhoh = (vol_c_rhoh./C_rhoh + Beta_rhoh)/sigma_b;

dC_dRho = (C_rhoh - C_b)/(rho_h - rho_b);
dB_dRho = (Beta_rhoh - Beta_b)/(rho_h - rho_b);
dK_dRho = (K_rhoh - K_b)/(rho_h - rho_b);
dEps_dRho = (epsilon_C_rhoh - epsilon_C_b)/(rho_h - rho_b);


%% INTERPOLANTS

dC_dE_FUN = griddedInterpolant(Y_mat,PHI_mat,dC_dE);
dB_dE_FUN = griddedInterpolant(Y_mat,PHI_mat,dB_dE);
dK_dE_FUN = griddedInterpolant(Y_mat,PHI_mat,dK_dE);
dEps_dE_FUN = griddedInterpolant(Y_mat,PHI_mat,dEps_dE);

dC_dL_FUN = griddedInterpolant(Y_mat,PHI_mat,dC_dL);
dB_dL_FUN = griddedInterpolant(Y_mat,PHI_mat,dB_dL);
dK_dL_FUN = griddedInterpolant(Y_mat,PHI_mat,dK_dL);
dEps_dL_FUN = griddedInterpolant(Y_mat,PHI_mat,dEps_dL);

dC_dR_FUN = griddedInterpolant(Y_mat,PHI_mat,dC_dR);
dB_dR_FUN = griddedInterpolant(Y_mat,PHI_mat,dB_dR);
dK_dR_FUN = griddedInterpolant(Y_mat,PHI_mat,dK_dR);
dEps_dR_FUN = griddedInterpolant(Y_mat,PHI_mat,dEps_dR);

dC_dRho_FUN = griddedInterpolant(Y_mat,PHI_mat,dC_dRho);
dB_dRho_FUN = griddedInterpolant(Y_mat,PHI_mat,dB_dRho);
dK_dRho_FUN = griddedInterpolant(Y_mat,PHI_mat,dK_dRho);
dEps_dRho_FUN = griddedInterpolant(Y_mat,PHI_mat,dEps_dRho);

%% PLOTS


paperposition = [0 0 8 4.5]*2/3;
definition = '-r500';
axis_font = 6;
label_font = 9;
legend_font = 7;

edgecolor = 'k';

linewidth = 1*2/3;

color1 = 'b';
color2 = 'r';
color3 = [0.4940 0.1840 0.5560];

marker1 = 'o';
marker2 = '^';
marker3 = '*';

marker_size = 5.5*2/3;

marker_spacing = [1:20:n_grid];



z1 = 0;
z2 = -0.25;
z3 = -0.5;


z1_lab = char(strcat('$$y$$ = ',{' '}, string(-z1)));
z2_lab = char(strcat('$$y$$ = ',{' '}, string(-z2)));
z3_lab = char(strcat('$$y$$ = ',{' '}, string(-z3)));


y1 = z1.*(1-PHI);
y2 = z2.*(1-PHI);
y3 = z3.*(1-PHI);

%% Eta
% C
VAR1 = dC_dE_FUN(y1,PHI);
VAR2 = dC_dE_FUN(y2,PHI);
VAR3 = dC_dE_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$dc/d\eta$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','southwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dC_dE_3line.png'],'-r500');
hold off; clear h; close;

% K
VAR1 = dK_dE_FUN(y1,PHI);
VAR2 = dK_dE_FUN(y2,PHI);
VAR3 = dK_dE_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$dk/d\eta$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','northwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dK_dE_3line.png'],'-r500');
hold off; clear h; close;


% Beta
VAR1 = dB_dE_FUN(y1,PHI);
VAR2 = dB_dE_FUN(y2,PHI);
VAR3 = dB_dE_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$d\hat\beta/d\eta$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','northwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dB_dE_3line.png'],'-r500');
hold off; clear h; close;

% Epsilon
VAR1 = dEps_dE_FUN(y1,PHI);
VAR2 = dEps_dE_FUN(y2,PHI);
VAR3 = dEps_dE_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$d\epsilon_C/d\eta$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','northwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dEps_dE_3line.png'],'-r500');
hold off; clear h; close;


%% Lambda
% C
VAR1 = dC_dL_FUN(y1,PHI);
VAR2 = dC_dL_FUN(y2,PHI);
VAR3 = dC_dL_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$dc/d\lambda$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','northwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dC_dL_3line.png'],'-r500');
hold off; clear h; close;

% K
VAR1 = dK_dL_FUN(y1,PHI);
VAR2 = dK_dL_FUN(y2,PHI);
VAR3 = dK_dL_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$dk/d\lambda$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','southwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dK_dL_3line.png'],'-r500');
hold off; clear h; close;


% Beta
VAR1 = dB_dL_FUN(y1,PHI);
VAR2 = dB_dL_FUN(y2,PHI);
VAR3 = dB_dL_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$d\hat\beta/d\lambda$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','southwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dB_dL_3line.png'],'-r500');
hold off; clear h; close;

% Epsilon
VAR1 = dEps_dL_FUN(y1,PHI);
VAR2 = dEps_dL_FUN(y2,PHI);
VAR3 = dEps_dL_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$d\epsilon_C/d\lambda$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','southwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dEps_dL_3line.png'],'-r500');
hold off; clear h; close;

%% R
% C
VAR1 = dC_dR_FUN(y1,PHI);
VAR2 = dC_dR_FUN(y2,PHI);
VAR3 = dC_dR_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$dc/dr$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','southwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dC_dR_3line.png'],'-r500');
hold off; clear h; close;

% K
VAR1 = dK_dR_FUN(y1,PHI);
VAR2 = dK_dR_FUN(y2,PHI);
VAR3 = dK_dR_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$dk/dr$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','southwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dK_dR_3line.png'],'-r500');
hold off; clear h; close;


% Beta
VAR1 = dB_dR_FUN(y1,PHI);
VAR2 = dB_dR_FUN(y2,PHI);
VAR3 = dB_dR_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$d\hat\beta/dr$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','northwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dB_dR_3line.png'],'-r500');
hold off; clear h; close;

% Epsilon
VAR1 = dEps_dR_FUN(y1,PHI);
VAR2 = dEps_dR_FUN(y2,PHI);
VAR3 = dEps_dR_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$d\epsilon_C/dr$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','northwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dEps_dR_3line.png'],'-r500');
hold off; clear h; close;


%% Rho
% C
VAR1 = dC_dRho_FUN(y1,PHI);
VAR2 = dC_dRho_FUN(y2,PHI);
VAR3 = dC_dRho_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$dc/d\rho$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','northwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dC_dRho_3line.png'],'-r500');
hold off; clear h; close;

% K
VAR1 = dK_dRho_FUN(y1,PHI);
VAR2 = dK_dRho_FUN(y2,PHI);
VAR3 = dK_dRho_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$dk/d\rho$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','southwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dK_dRho_3line.png'],'-r500');
hold off; clear h; close;


% Beta
VAR1 = dB_dRho_FUN(y1,PHI);
VAR2 = dB_dRho_FUN(y2,PHI);
VAR3 = dB_dRho_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin-0.5 ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$d\hat\beta/d\rho$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','southwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dB_dRho_3line.png'],'-r500');
hold off; clear h; close;

% Epsilon
VAR1 = dEps_dRho_FUN(y1,PHI);
VAR2 = dEps_dRho_FUN(y2,PHI);
VAR3 = dEps_dRho_FUN(y3,PHI);
ymin = min([nanmin(VAR1), nanmin(VAR2), nanmin(VAR3)]);
ymax = max([nanmax(VAR1), nanmax(VAR2), nanmax(VAR3)]);
figure;
hold on
l1 = plot(PHI, VAR1, 'Color', color1, 'LineWidth', linewidth, 'Marker', marker1, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l2 = plot(PHI, VAR2, 'Color', color2, 'LineWidth', linewidth, 'Marker', marker2, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
l3 = plot(PHI, VAR3, 'Color', color3, 'LineWidth', linewidth, 'Marker', marker3, 'MarkerSize', marker_size, 'MarkerIndices', marker_spacing);
ylim([ymin ymax])
set(gca,'FontSize',axis_font);
Xla=xlabel('Beliefs ($$\phi$$)','Interpreter','Latex','FontSize',label_font);
set(Xla, 'Units', 'Normalized');
Yla=ylabel('$$d\epsilon_C/d\rho$$','Interpreter','Latex','FontSize',label_font);
set(Yla, 'Units', 'Normalized');
set(gcf,'PaperUnits','centimeters')
set(gcf, 'PaperPosition', paperposition);
[hh,icons,plots,txt] = legend([l1 l2 l3],{z1_lab,z2_lab, z3_lab},'Interpreter', 'latex', 'Location','southwest', 'FontSize', legend_font);
icons(4).LineStyle = 'none';
icons(1).Position = [icons(1).Position(1)*1/4*1.2 icons(1).Position(2) icons(1).Position(3)];
icons(5).XData = icons(5).XData/4;
icons(2).Position = [icons(2).Position(1)*1/4*1.2 icons(2).Position(2) icons(2).Position(3)];
icons(6).LineStyle = 'none';
icons(7).XData = icons(7).XData/4;
icons(3).Position = [icons(3).Position(1)*1/4*1.2 icons(3).Position(2) icons(3).Position(3)];
icons(8).LineStyle = 'none';
icons(9).XData = icons(9).XData/4;
legend boxoff
print('-dpng',gcf,['dEps_dRho_3line.png'],'-r500');
hold off; clear h; close;